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Abstract 

Geometric  design,  including  the  internal  structure  and  external  shape,  considerably  affect  the  thermal,  fluid,  and  electrochemical  characteristics 
of  a  polymer  electrolyte  membrane  (PEM)  fuel  cell,  which  determine  the  polarization  curves  as  well  as  the  thermal  and  power  inertias.  Shape 
optimization  is  a  natural  alternative  to  improve  the  fuel  cell  performance  and  make  fuel  cells  more  attractive  for  power  generation.  This  paper 
investigates  the  internal  and  external  structure  effects  on  the  fuel  cell  steady  and  transient  operation  with  consideration  of  stoichiometric  ratios, 
pumping  power,  and  working  temperature  limits.  The  maximal  steady  state  net  power  output  and  the  fuel  cell  start-up  time  under  a  step-changed 
current  load  characterize  the  fuel  cell  steady  and  transient  performance  respectively.  The  one-dimensional  PEM  fuel  cell  (PEMFC)  thermal  model 
introduced  in  a  previous  work  [J.V.C.  Vargas,  J.C.  Ordonez,  A.  Bejan,  Constructal  flow  structure  for  a  PEM  fuel  cell,  Int.  J.  Heat  Mass  Transfer 
47  (2004)  4177-4193]  is  amended  to  simulate  the  fuel  cell  transient  start-up  process.  The  shape  optimization  consists  of  the  internal  and  external 
PEMFC  structure  optimization.  The  internal  optimization  focuses  on  the  optimal  allocation  of  fuel  cell  compartment  thicknesses.  The  external 
optimization  process  seeks  the  PEM  fuel  cell  optimal  external  aspect  ratios.  These  two  levels  of  optimizations  pursue  the  optimal  geometric  design 
with  quick  response  to  the  step  loads  and  large  power  densities.  Appropriate  dimensionless  groups  are  identified  and  the  numerical  results  are 
presented  in  dimensionless  charts  for  general  engineering  design.  The  universality  of  the  general  optimal  shape  found  is  also  discussed. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Fuel  cells  are  expected  to  become  of  common  use  in  sta¬ 
tionary  power  generation  and  transportation  systems  [1-3].  The 
energy  industry  is  searching  for  alternative  power  generation 
modes  in  order  to  reduce  its  dependency  on  conventional  fossil 
fuel  and,  in  parallel,  for  means  to  improve  systems  and  processes 
efficiencies.  Fuel  cell  technology  also  benefits  the  environment 
by  cutting  down  a  large  amount  of  green  house  gas  (CO2)  and 
other  poisonous  exhausts  (NOx,  SOx),  through  the  introduc¬ 
tion  of  the  fuel  cells  into  transportation  and  power  generation 
systems.  Furthermore,  submarines  and  aircrafts  with  fuel  cell 
power  supply  system  could  dramatically  decrease  cruising  noise 
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to  comply  with  military  requirements.  PEMFC  is  also  widely 
investigated  because  of  its  high  efficiency  and  the  potential  for 
utilization  in  vehicular  and  portable  applications.  Fuel  cell  tech¬ 
nology  has  been  improving  in  several  aspects  such  as  materials, 
thermodynamic  management  and  flow  structure.  In  this  paper, 
we  focus  on  the  flow  structure  optimization,  which  has  the  poten¬ 
tial  to  improve  the  fuel  cell  performance,  involving  steady  state 
and  transient  operation. 

The  flow  structure  optima  are  the  results  of  numerical  simu¬ 
lation  of  mathematical  models  that  express  the  physics  of  fuel 
cell  operations.  Kumar  and  Reddy  [4]  present  a  mathematic 
model  applying  Navier-Stokes  equations  to  predict  the  effects 
of  shape  and  size  of  flow  channel  on  the  fuel  cell  performance. 
Yi  and  Nguyen  [5]  developed  a  two-dimensional  cathode  model 
to  investigate  the  performance  of  interdigitated  fuel  cells.  Zhou 
and  Liu  [6]  presented  a  three-dimensional  steady  state  model  for 
PEM  fuel  cells,  Hu  et  al.  [7,8]  introduced  a  three-dimensional 
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Nomenclature 

A  a 

variation  of  quantity  a  between  two  conditions 

A 

area  (m2) 

A 

dimensionless  area 

As 

fuel  cell  cross-sectional  area  (m2) 

Aw 

the  wall  area  (m2) 

^wet 

the  wetted  surface  area  (m2) 

B 

dimensionless  constant  in  Eq.  (14) 

cp 

specific  heat  at  constant  pressure  (kJkg-1  K-1) 

Oy 

specific  heat  at  constant  volume  (kJkg-1  K-1) 

cv 

control  volume 

D 

Knudsen  diffusion  coefficient  (m2  s-1) 

Dh 

gas  channel  hydraulic  diameter  (m) 

E 

energy  flow  rate  (W) 

E 

dimensionless  energy  flow  rate 

f 

friction  factor 

F 

faraday  constant  (96,500  Cequiv.-1) 

G 

molar  Gibbs  free  energy  change  (kJkmol-1) 

G 

dimensionless  Gibbs  free  energy  change 

h 

heat  transfer  coefficient  (W  m-2  K-1) 

h 

dimensionless  heat  transfer  coefficient 

Hi(Ti ) 

molar  enthalpy  of  formation  at  a  temperature  7) 
of  reactants  and  products  (kJkmol-1) 

Him 

dimensionless  molar  enthalpy  of  formation  at  a 
dimensionless  temperature  0*  of  reactants  and 
products 

lo,c 

exchange  current  densities  (Am-2) 

/lima,  /iim,c  limiting  current  densities  (Am  2) 

/ 

current  (A) 

I 

dimensionless  current 

j 

mass  flux  (kg  s-1  m-2) 

k 

thermal  conductivity  (W  m-1  K-1) 

k 

dimensionless  thermal  conductivity 

K 

permeability  (m2) 

L 

length  (m) 

Lchi  Lt 

gas  channels  internal  dimensions  as  shown  in 
Fig.  1  (m) 

Lx,  Ly,  Lz  fuel  cell  length,  width  and  height,  respectively 
(m) 

m 

mass  flow  rate  (kg  s-1) 

M 

molecular  weight  (kgkmol-1) 

n 

equivalent  electron  per  mole  of  reactant 
(equiv.  mol-1) 

h 

molar  flow  rate  (kmol  s-1) 

^ch 

number  of  parallel  ducts  in  gas  channel 

N 

dimensionless  global  wall  heat  transfer  coefficient 

P 

pressure  (N  m-2) 

P 

dimensionless  pressure 

PEMFC 

polymer  electrolyte  membrane  fuel  cell 

Pr 

Prandtl  number  (p xp  K-1) 

a 

tortuosity 

Q 

reaction  quotient 

Q 

heat  transfer  rate  (W) 

Q 

dimensionless  heat  transfer  rate 

r 

pore  radius  (m) 

R 

ideal  gas  constant  (kJkg-1  K-1) 

R 

universal  gas  constant  8.314  (kJkmol-1  K-1) 

Re  h 

Reynolds  number  based  on  (Dh,  uD^pl/i) 

s 

dimensionless  conversion  factor,  Eq.  (37) 

t 

time  (s) 

T 

temperature  (K) 

u 

mean  velocity  (ms-1) 

u 

dimensionless  mean  velocity 

U 

global  wall  heat  transfer  coefficient  (W  m-2  K- 1 ) 

V 

voltage  (V) 

Vt 

total  volume  (m3) 

V 

dimensionless  voltage 

Vt 

dimensionless  total  volume 

W 

dimensionless  fuel  cell  electrical  power 

Wne, 

dimensionless  fuel  cell  net  power,  Eq.  (33) 

^net,s 

dimensionless  fuel  cell  steady  state  net  power  out¬ 
put 

tip 

dimensionless  required  pumping  power,  Eqs.  (33) 
and  (36) 

[] 

molar  concentration  of  a  substance  (moll-1) 

Greek  symbols 

C^a, 

anode  and  cathode  charge  transfer  coefficients 

/? 

electrical  resistance  (Q) 

y 

specific  heat  ratio 

8 

gas  channel  aspect  ratio 

? 

stoichiometric  ratio 

*7a»  0c 

anode  and  cathode  charge  transfer  overpotentials 
(V) 

0d,a->  0d,c 

anode  and  cathode  mass  diffusion  overpotentials 
(V) 

0 a,  be 

dimensionless  anode  and  cathode  charge  transfer 
overpotentials 

i)d,a?  i?d,< 

:  dimensionless  anode  and  cathode  mass  diffusion 
overpotentials 

7  ohm 

dimensionless  fuel  cell  total  ohmic  potential  loss 

0 

dimensionless  temperature 

X 

ionomer  water  content 

li 

viscosity  (Pa  s) 

V 

reaction  coefficients 

% 

dimensionless  length 

0 

porosity 

P 

density  (kg  m-3) 

P 

dimensionless  density 

O 

electrical  conductivity  (£2-1  m-1) 

X 

dimensionless  time 

iff 

dimensionless  mass  flow  rate 

Subscripts 

a 

anode 

(aq) 

aqueous  solution 

c 

cathode 

cond 

heat  conduction 
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conv  heat  convection 

CV  control  volume 

e  reversible 

f  fuel 

(g)  gaseous  phase 

h  hydraulic 

H+  hydrogen  cation 

H2  hydrogen 

H2O  water 

in  control  volume  inlet 

j  control  volume  index 

k  control  volume  index 

(1)  liquid  phase 

m  one-way  maximum 

mm  two-way  maximum 

ohm  ohmic 

opt  optima 

out  control  volume  outlet 

ox  oxidant 

O2  oxygen 

p  polymer  electrolyte  membrane 

ref  reference 

sa  anode  solid  side 

sc  cathode  solid  side 

su  start-up 

w  wall 

wet  wetted  surface 

0  initial  condition 

1,  . . 7  control  volume  index 
12  interaction  between  CV  1  and  CV2 

23  interaction  between  CV2  and  CV3 

34  interaction  between  CV3  and  CV4 

45  interaction  between  CV4  and  CV5 

56  interaction  between  CV5  and  CV6 

67  interaction  between  CV6  and  CV7 

00  ambient 

Superscript 

0  standard  conditions  [gases  at  1  atm,  25  °C,  species 

in  solution  at  1 M,  where  M  is  the  molar¬ 
ity  =  (moles  solute)/(liters  solution)] 


two  phase  flow  mathematical  model  with  consideration  of  exist¬ 
ing  multi-phase  water  in  the  flow  fields.  Those  models  are  expen¬ 
sive  to  solve  numerically,  due  to  the  partial  difference  equations, 
and  to  the  large  number  of  flow  configurations  to  be  tested  in  an 
optimization  procedure.  At  the  PEMFC  stack  level,  Amphlett 
et  al.  [9]  introduced  a  transient  mathematical  model  for  predict¬ 
ing  the  response  and  performance  of  the  system  for  applications 
where  the  operating  conditions  change  with  time,  considering 
heat  losses,  changes  in  stack  temperature,  reactant  gas  concen¬ 
trations,  and  other  internal  phenomena,  but  no  consideration  was 
given  to  the  pressure  drops  experienced  by  the  fuel  and  oxidant 
in  the  feeding  headers  and  gas  channels.  Experimental  results 


were  presented  by  Hamelin  et  al.  [10]  for  the  time  response  of 
a  stationary  PEMFC  stack  under  fast  load  commutations,  and 
compared  to  the  Amphlett  et  al.  [9]  model  with  good  agree¬ 
ment.  Several  other  experimental  studies  have  been  conducted 
to  access  the  dynamic  behavior  of  PEMFC  stacks  in  search  for 
performance  improvement  in  the  transient  mode  [11-15].  Recent 
studies  presented  PEMFC  stack  steady  state  mathematical  mod¬ 
els  to  predict  performance  [16-18]:  one  study  [18]  presented  a 
genetic  algorithm  technique  for  finding  the  best  configuration  of 
the  stack  in  terms  of  number  of  cells  and  cell  surface  area,  but  did 
not  consider  headers  and  gas  channels  pumping  power  losses. 
No  transient  model  was  found  in  the  literature  that  addresses  the 
spatial  temperature  and  pressure  gradients  in  a  PEMFC  stack, 
pressure  drops  in  the  headers  and  single  cells  gas  channels  and 
their  effect  on  performance,  with  the  exception  of  a  recent  study 
published  by  Vargas  et  al.  [19]. 

As  a  result,  for  the  fuel  cell  structure  optimization,  it  is  desir¬ 
able  that  the  model  is  as  simple  as  possible,  in  order  to  require 
low  computational  time  to  provide  solutions  for  each  tested  con¬ 
figuration,  thus  allowing  for  an  effective  optimization  procedure. 
With  such  objectives  in  mind,  Vargas  and  Bejan  [20],  and  Vargas 
et  al.  [21]  presented  a  lumped  control  volume  model  for  alka¬ 
line  and  PEM  fuel  cells,  respectively,  which  calculated  the  fuel 
cell  thermal-fluid  field,  mass  transport  mechanism,  electrochem¬ 
ical  properties,  and  power  consumption  and  generation.  The 
optimal  fuel  cell  geometry,  including  internal  layer  thickness 
allocation  and  external  aspect  ratios,  were  also  obtained  in  those 
studies. 

In  this  paper,  we  amend  the  model  introduced  by  Vargas  et 
al.  [21]  and  extend  the  flow  structure  optimization  procedure  to 
multi-objective  PEMFC  optimization  involving  the  steady  state 
and  transient  behaviors.  The  spatial  temperature  and  pressure 
gradients  in  a  single  PEMFC,  pressure  drops  in  the  cell  gas 
channels,  and  their  effect  on  performance  are  investigated.  The 
multi-objective  optima  according  to  constructal  theory  [22]  are 
also  discussed  in  this  paper. 


2.  Thermodynamic  model 

A  thermo-electrochemical  model  for  a  single  PEM  fuel  cell 
operating  at  steady  state  has  been  introduced  in  a  previous  study 
by  Vargas  et  al.  [21].  In  this  study,  that  mathematical  model 
is  amended  to  study  the  PEM  fuel  cell  transient  performance. 
According  to  Vargas  et  al.  [21],  the  fuel  cell  is  divided  into 
seven  control  volumes  that  are  shown  in  Fig.  1.  These  control 
volume  (CV)  interact  energetically  with  one  another.  The  fuel 
cell  also  interacts  with  adjacent  fuel  cells  in  a  stack,  and/or 
with  the  ambient.  Additionally,  two  bipolar  plates  (intercon¬ 
nects)  have  the  function  of  allowing  the  electrons  produced 
by  the  electrochemical  oxidation  reaction  at  the  anode  to  flow 
to  the  external  circuit  or  to  an  adjacent  cell.  The  control  vol¬ 
umes  are  the  fuel  channel  (CV1),  the  anode  diffusion-backing 
layer  (CV2),  the  anode  reaction  layer  (CV3),  the  polymer  elec¬ 
trolyte  membrane  (CV4),  the  cathode  reaction  layer  (CV5),  the 
cathode  diffusion  backing  layer  (CV6)  and  the  oxidant  channel 
(CV7). 
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Fig.  1.  PEM  fuel  cell  internal  structure  and  external  shape. 


The  model  consists  of  the  conservation  equations  for  each 
control  volume,  and  the  equations  accounting  for  electrochem¬ 
ical  reactions,  where  reactions  are  present.  The  reversible  elec¬ 
trical  potential  and  power  of  the  fuel  cell  are  then  computed 
(based  on  the  reactions)  as  functions  of  the  temperature  and 
pressure  fields  determined  by  the  model.  The  actual  electri¬ 
cal  potential  and  power  of  the  fuel  cell  are  obtained  by  sub¬ 
tracting  from  the  reversible  potential  the  losses  due  to  surface 
overpotentials  (poor  electrocatalysis),  slow  diffusion  and  all 
internal  ohmic  losses  through  the  cell  (resistance  of  individ¬ 
ual  cell  components,  including  electrolyte  layer,  interconnects 
and  any  other  cell  components  through  which  electrons  flow). 
These  are  functions  of  the  total  cell  current  (7),  which  is  directly 
related  to  the  external  load  (or  the  cell  voltage).  In  sum,  the 
total  cell  current  is  considered  an  independent  variable  in  this 
study. 

2.7.  Dimensionless  variables 

To  make  the  results  general,  this  model  and  its  simulation 
results  are  formulated  and  reported  in  a  dimensionless  way.  The 
dimensionless  variables  are  defined  based  on  the  geometric  and 
operating  conditions  of  the  PEM  fuel  cell. 


The  dimensionless  length  §  and  dimensionless  area  A  are 
defined  as 


L 


V 


1/3  ’ 


T 


A  = 


A 


Vr 


2/3 


T 


where  Vr  =  LxLyLz,  is  the  total  fuel  cell  volume. 

Pressures  and  temperatures  are  referred  to  ambient  condi¬ 
tions:  P = p/poo  and  0  =  TIT0 0 .  The  dimensionless  mass  flow  rates 
are  defined  as  follows: 


where  the  subscript  i  indicates  the  reactant  or  product  that  flows 
through  the  fuel  cell;  and  mref  is  a  specified  reference  mass  flow 
rate  given  in  Table  1.  The  reference  time  fref  depends  on  the 
reference  mass  flow  rate  mre f,  and  therefore,  the  dimensionless 
time  r  and  reference  time  tTQ f  are  expressed  as  follows: 


^ref 


Poo  Tt 


tffMref 


where  the  subscript  “f”  stands  for  fuel.  In  this  paper,  for  sim¬ 
plicity,  we  assume  that  the  fuel  stream  is  pure  hydrogen,  and  the 
oxidant  is  pure  oxygen. 
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Table  1 

Properties  and  constant  values  used  in  the  simulations 
5  =  0.156 

Cp,f=  14.95  kJ  kg-1  K-1 
Cpox  =  0.91875  kj  kg-1  K_1 
c„,f  =  10.8  kj  kg-1  K-1 
cVfOX  =  0.659375  kj  kg"1  K_1 
Cv, sa  =  c„,sc  =  133  J  kg-1  K-1 
CSoi  =  4180  Jkg-1  K-1 
F = 96,500  C  mol-1 
<o,a  =  <’o,c  =  10  A  m-2 
Act  =  1  A 

kf  =  0.2  W  m1  K-1 
kox  =0.033  W  m-1  K"1 
Ap  =  0.21  W  m_1  K"1 
As,a  =  As,c  =  7 1 .6  W  m_1  K-1 
K2  =  K6  =  4  x  10“ 14  m2 
K3  =  K,=4x  10“16m2 
mre  f  =  10-4kgs_1 
Pi'  =  Poo  =  o.  1  MPa 
Pox  =  0.12  MPa 
q  =  1.5 

5f  =  4.157  kJ  kg-1  K-1 
Rox  =  0.2598  kJ  kg-1  K-1 
Tf  =  Tox  =  Too  =298.15  K 
Uwi  =  50  Wm“2  K,  i  =  l,7 
Vref=lV 

Vf  =  2.25  x  10-5  m3 

VT,ref=10-5m3 

O' a  =  Q'c  =  0.5 

0  2=06=0.4 

03  =05  =0.2 

pi  =  10“5  Pas 

/X7  =2.4  x  1(T5  Pas 

o'!,  07  -  1.388  x  106  m_1 

(72,  cr 5  =  4000  m_1 


The  physics  of  the  fuel  cell  is  described  by  taking  into 
account  the  mass  conservation  and  the  first  law  of  thermody¬ 
namics  at  each  CV,  and  the  electrochemical  reactions  at  CV3  and 
CV5. 


2.2.  Mass  balance 


The  hydrogen  mass  flow  rate  required  by  the  working  current 
(7),  which  indicates  the  external  load,  is  given  by 

7 

^h2  =  5h2  Mh2  =  — A7h2  (7) 

nF 

where  Ah2  is  the  molar  flow  rate  for  hydrogen;  Mh2  the  molec¬ 
ular  weight  of  hydrogen;  n  the  equivalent  electron  per  mole  of 
reactant,  equiv.  mol-1  (n  =  2);  and  F  is  the  Faraday  constant, 
96,500  Cequiv.-1. 

Similarly,  the  oxygen  mass  flow  rate  needed  in  the  PEM  fuel 
cell  electrochemical  reaction  is  evaluated  as  follows: 


1 

mo2  =  2m2Mo2 


The  stoichiometric  ratios,  denoted  as  f  i  on  the  fuel  side  and 
£7  on  the  oxidant  side,  are  defined  as  the  provided  reactant 
divided  by  the  reactant  needed  for  the  electrochemical  reac¬ 
tion  of  interest  and  greater  than  1.  The  stoichiometric  ratios 
are  the  design  constraints  of  this  study.  Therefore,  considering 
the  stoichiometric  ratios,  the  dimensionless  inlet  reactant  mass 
flow  rates  are  Vt  =  ?i^h2  =  £i^H2Mref  and  x/r ox  =  ?7^02  = 
fimo2/mVQf,  where  mH2and  fho2  are  obtained  from  Eqs.  (7) 
and  (8). 


2.3.  Energy  conservation 


Other  variables  are  the  nondimensional  global  wall  heat  trans¬ 
fer  coefficient  N,  heat  transfer  coefficient  h,  thermal  conductivity 
k,  density  p  and  specific  heat  ratio  y,  defined  as  follows: 


N  = 


P  = 


UyV- F 

FlrefC  pf 
pRfToo 
Pf 


hv 2^3 
h  =  — T- 


kv 

~k=JWj_ 


y  = 


cp,  f 


c 


(4) 


V 


where  cpy  is  the  fuel  specific  heat  at  constant  pressure. 

The  reference  energy  flow  rate  is  defined  by  Eref  = 
threfCpfToo,  which  is  used  as  a  reference  scale  for  power  related 
variables: 


(A E,  AECW,  Q,  AH,  AG)  = 


(A E,  AECw,  Q,  AH,  AG) 

E  ref 


(5) 


The  dimensionless  heat  transfer  rates  are  given  by 

—kAAO 


Qcoxw  —  hAAO, 

72 

Q  ohm  = 


2cond  — 


FlrefCpfToQ 

where  ft  is  the  electrical  resistance. 


0W  =  NAwA0, 


(6) 


The  thermal  model  consists  of  coupled  ordinary  differ¬ 
ential  equations  that  are  obtained  from  the  energy  conser¬ 
vation  principle  applied  to  each  control  volume,  as  shown 
in  Fig.  1.  The  energy  conservation  in  CV1  is  expressed  as 
follows: 


Q w  1  T  ^l)  +  012  +  f2ohml  — 


1  dE’cvi 


EfQf  d  t 


(9) 


where  Qw\  =  AfiAwi(l  —  0\),  stands  for  the  heat  transfer  rate 
from  the  ambient  to  the  CV  1  through  the  wall.  The  wall  surface 
connecting  the  CV)  (j=  1,  7)  and  the  ambient  (AWJ)  is  calcu¬ 
lated  by  AWj  =  +  2 (£y  +  t;z)i;j.  The  second  term  in  the  left 

hand  side  of  Eq.  (9)  is  the  dimensionless  fuel  enthalpy  decrease, 
which  is  obtained  by  assuming  that  the  reactant  leaves  the 
fuel  cell  at  the  CV  1  temperature.  Qn  =  h\As(l  —  </>2)($2  —  0\) 
denotes  the  heat  transfer  rate  from  CV2  to  CV1.  The  fuel 
cell  cross-sectional  area  is  As  =  The  heat  transfer  coef¬ 
ficients  hj  (j=  1,  7)  in  the  gas  channels  are  estimated  accord¬ 
ing  to  the  flow  regime.  For  the  laminar  flow  (7feh<2300) 
[23]: 


hjDh,j 


=  7.541(1  -  2.6108 j  +  4.970S2  -5.1195 


* J 


j 


+2.702^-0.548^-),  j  =  1,7 


(10a) 
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For  the  turbulent  flow  (2300  <  Re h  <  2  x  104)  [24] : 


(fj/2)(Rehj  -  1000) /V; 

1  +  l2.1(fj/2)l/2(Pr2f  -  1)’ 


(10b) 


where  Sj  =  L^/Lj,  for  Lch  <  Lj  and  Sj-Lj/Lch,  for  Lch  >  Lj ; 
Dh,/  —  2Lc]]L//( Lc\]  +  Fy),  Pv  is  the  gas  Prandtl  number;  and 

/  is  the  friction  factor  given  by  fj  =  0.079  Reh  V4  when 

2300  <  Re^j  <  2  x  104  [24]. 

The  Reynolds  number  is  a  function  of  the  gas  flow  velocity 
Uj.  According  to  Vargas  et  al.  [21],  these  velocities  are  calculated 
by 


u  i 


PfLchL\Hch 


(Ha) 


W7 


RoxTjPloX  f  , 
2poxLc\\LlWch  \ 


(11b) 


Finally,  the  ohmic  heating  rate  <20hmi  =  /2 ^i/(^ref^,f^oo) 
with  the  previous  dimensionless  energy  fluxes  contribute  to  the 
increase  of  the  interconnect  internal  energy  as  follows: 


dOi 

dr 


Yf  01,0 
Pi 


[0wl  +  VT(0f  “  #l)  +  012  +  0ohml] 


(12) 


where  j2  —  7hH2/A3?wet>  and  j3  —  ^o2/^5,wet-^l3,wet  andA5  ?wet 
are  the  wetted  surface  areas  in  the  porous  catalyst  layers,  defined 

as  Awetj  =  4<fijLjKj  As,  j  =  3,  5,  where  K  is  the  perme¬ 
ability  and  0  is  the  porosity.  Because  the  channel  layer  is  struc¬ 
turally  connected  to  CV2  and  CV6,  we  assume  that  P2,m  =  Pf 
and  P6,in  =  Fox-  The  average  pressures  in  CV2  and  CV6  are  esti¬ 
mated  as  follows: 

Pj  —  2  (T* in  +  P 7, out)?  j  =  2,  6  (16) 

In  the  reaction  layer  (CV3),  there  are  a  few  consider¬ 
ations  that  must  be  taken  into  account  in  the  energy  bal¬ 
ance.  An  oxidation  reaction,  Ffyg)  — >  2H(aq)+  +  2e_,  takes  place 
in  this  layer.  The  dimensionless  molar  enthalpy  change  due 

to  the  reaction  is  given  by  -A  H3  =  E reactants  I '  A/7/)  I  ~ 
^products [y7^7'(T/)].  Here,  Vj  is  the  reaction  coefficient.  Hj  is 
the  dimensionless  molar  enthalpy  of  formation  of  component 
j  [28,29].  The  dimensionless  reversible  electrical  power,  WQ 3, 
is  expressed  in  terms  of  the  dimensionless  Gibbs  free  energy 
change,  i.e.,  We3  =  —  AG 3,  which  is  also  part  of  the  CV3’s 
energy  balance.  The  Gibbs  free  energy  change,  AG3,  is  a  func¬ 
tion  of  the  temperature  and  pressure  as  follows: 

A Gj  =  A G°j  +  RTj  In  Qj,  j  =  3,  5  (17) 


where  6 i,o  is  the  CVl’s  initial  dimensionless  temperature,  and 
Yf  =  cp  f/ cv  {  is  the  fuel  gas  specific  heat  ratio.  The  average 
fuel  pressure  (pf  in  CV1)  and  oxidant  pressure  (pox  in  CV7)  are 
assumed  known  and  constants,  shown  in  Table  1,  during  the  fuel 
cell  operation. 

The  anode  backing  diffusion  layer  (CV2)  is  porous,  which 
consists  of  a  solid  side  and  a  fluid  side.  The  mass  of  fluid  in  CV2 
is  negligible  compared  to  the  mass  of  solid.  Therefore,  the  CV2 
internal  energy  differential  is  d£cv2  =  ^cv2G,sad72.  The  net 
heat  transfer  rate  through  CV2  is  Q2  =  —  012  +  0W2  +  023  + 
2ohm2,  where  Q2 3  =  -£sa(  1  -  02)AS(6>2  -  03)/[(£2  +  ?3)/2]. 
The  anode  solid  side  dimensionless  thermal  conductivity  ksa 
was  defined  in  Eq.  (4).  The  energy  balance  in  CV2  leads  to  the 
CV2’s  dimensionless  temperature  expression  as  follows: 


d02 

dr 


Psa(l 


Ysa 

~  02)§2 SySz 


[0H2(£l 


@2)  +  02] 


(13) 


where  psais  the  dimensionless  density  of  the  anode  solid  side, 
and  ysa  is  the  effective  specific  heat  ratio  of  the  anode  solid  side. 

The  flow  in  the  electrodes  is  modeled  as  Knudsen  flow  [25]. 
The  fuel  and  oxidant  mass  fluxes  are  given  by 


jj  —  [£>(Pout  Pin)/F]y,  j  —  2,  6  (14) 

where  D  =  B{r[8RT/(irM)l/2(pq}  is  the  Knudsen  diffusion 
coefficient;  p  the  density;  R  the  universal  gas  constant;  0  the 
porosity;  q  the  tortuosity  [26,27];  B  is  a  correction  coefficient. 
By  using  Eq.  (14)  and  the  ideal  gas  model  for  H2  and  O2,  we  find 
the  pressures  of  the  hydrogen  and  oxygen  that  enter  the  catalyst 
layers  as  follows: 


j j  RkTooL  ft  j 

7)  j  Poo 


j  =  2,  6;  k  =  f,  ox 


(15) 


where  AG3  is  the  standard  Gibbs  free  energy  change 
(kJkmol-1);  Q  is  the  reaction  quotient  and  is  formulated  by 

£>3  =  [/7(aq)+]  /  Pu2 ;  [/7(aq)+]  is  the  molar  concentration  of  the 
acid  solution  (mol/1);  and  pu2  =  p2,out  is  the  partial  pressure 
of  H2  at  the  CV2  outlet.  Recall  that  the  pure  liquids  or  solids  do 
not  appear  in  the  calculation  of  Q3 ,  neither  does  the  solvent  in 
a  dilute  solution,  [H(aq)+]  is  a  function  of  water  content  A,  i.e., 
[7/(aq)+]  ~  Ph2o / 7-4 47h2o  for  a  dilute  water  solution. 

The  dimensionless  net  heat  flux  flowing  in  CV3  is  given  by 
03  =  -  023  +  0w3  +  034  +  0ohm3  •  Since  the  heat  transfer 
rate  between  CV3  and  CV4  (the  polymer  electrolyte  membrane) 
is  dominated  by  conduction,  therefore,  this  dimensionless  heat 
flux  is  034  =  -2(1  -  03)(P3  -  P4)Asksakp/(§4ksa  +  §3kp). 

Considering  the  solid  compartment  dominates  the  electrode 
mass,  the  mass  and  energy  balance  analysis  in  CV3  leads  to  a 
dimensionless  temperature  expression  as  follows: 


dP3 

dr 


PsaQ 


Ysa 

-  03)§3 


Afo  +  AG3] 


(18) 


where  (A H3,  A G3)  =  hu2(AH3,  AG3)/(mrefCp^T00),  hu2  = 
mH2/AfH2,  and  other  parameters  are  similar  to  CV2  discussion. 

The  polymer  electrolyte  membrane  (CV4)  interacts  with 
CV3,  CV5  and  the  ambient.  In  the  cathode  reaction  layer 
(CV5),  there  occurs  the  reaction  expressed  by  502(g)  +  2e-  + 
2H(aq)+  — >  H20(i).  Considering  the  mass  conservation  in  CV4 
and  other  control  volumes  in  the  fuel  cell,  we  have  2Ph2  = 
/i H+  =  fiu+  =  2ho2.  In  conclusion,  ho2  =  hu2,  where  ho2  = 
2mo2/ Mo2.  Accordingly,  the  required  oxidant  mass  flow  rate 
is  rho2  =  mH2Mo2/2MH2.  The  dimensionless  net  heat  transfer 
in  CV4  is  obtained  from  Q4  =  -034  +  0w4  +  045  +  0ohm4 
and  045  =  2(1  05)(P4  ^5)Asksakp/(^4ksc  +  ^5^p)«  Then, 
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the  CV4  dimensionless  temperature  is  obtained  from 


d04 

dr 


PpU&Z 


04  +  H(63)H{m)+ 


H(04  )%q)+ 


(19) 


where  pp  is  the  dimensionless  density  of  membrane  solution,  and 
yp  is  the  effective  specific  heat  ratio  of  the  polymer  electrolyte 
membrane. 

The  analysis  in  the  cathode  reaction  layer  (CV5)  is  analogous 
to  what  we  previously  presented  in  the  anode  reaction  layer 
(CV3)  analysis.  The  CV5  dimensionless  temperature  is  obtained 
by 


d0s 

dr 


Psc(l 


Ysc 

-  05)§5 $y$z 


AH5  +  A  GS] 


(20) 


where  (A H5,  A G5)  =  ho2(AH5,  AG5)/(mrefcp^T00),  ho2  = 
mn?/2MH2,  fiH2o,out  =  fio2Psc  is  the  dimensionless  density  of 
the  cathode  solid  side,  and  ysc  is  the  effective  specific  heat  ratio 
of  the  cathode  solid  side. 

Similarly,  the  dimensionless  net  heat  transfer  rate  flow¬ 
ing  in  CV5  is  given  by  05  =  -  Q45  +  0w5  +  £>56  + 

2ohm5?  with  056  =  2kscAs(l  06005  06)/(§5  T  §6)- 

The  molar  enthalpy  change  (AH5  =  J2pmducts[vjHj(Tj)]- 

Sreactants [vjHj(Tj)])  indicates  the  reaction  heat  production 
(kJkmol-1);  Vj  are  the  reaction  coefficients.  14/(7))  is  the  for¬ 
mation  enthalpy  (kJkmoU1)  at  a  temperature  7)  of  compound 
j  in  reactants  and  products.  The  enthalpies  of  formation  are 
obtained  from  tabulated  values  [28,29]  at  for  02(g)  and 
T4  for  H(aq)+,  and  T$  for  H20(i)  at  1  atm,  respectively.  The 

change  in  the  Gibbs  free  energy,  AG5,  is  calculated  by  Eq.  (17). 

,  2  1/2  1 

The  CV5  reaction  quotient  is  05  =  {[77(aq)+]  p^}  ,  where 

P02  =  7*6,  out* 

The  mass  balance  for  CV6  yields  mo2  out  =  ^o2in  =  mo2  and 

hu2 o  =  ^H2o,out  =  ^H?o,in  =  fio2-  The  dimensionless  net  heat 
transfer  rate  in  CV6  results  from  06  =  —  056  +  0w6  +  067  + 

0ohm6?  067  =  ^7^s(l  ~  06007  “  0$),  h'j  =  hjVj  /mrQfCpj. 

The  dimensionless  temperature  for  CV6  is  given  by 


d06  _  Ysc 

Psc^Ty£z^6(l  06) 


X 


06  +  00,— (07 

CP,  f 


06)  +  ^&(05)h2O 


^(06)h2O 

(21) 


The  dimensionless  net  heat  transfer  rate  in  CV7  is  07  = 
—  067  +  0w7  +  0ohm6-  The  balances  for  mass  and  energy  in 
the  oxidant  channel  (CV7),  with  the  assumptions  of  non-mixing 
flow  and  the  assumption  that  the  space  is  filled  mainly  with  dry 

oxygen,  yields  mH2o  =  ^H2o,in  =  ^H2o,out  =  fio2+^H2o,  and 
the  dimensionless  CV7  temperature  as  follows: 


fi07  Xox  ^ox07,O 
dt  P ox^c^l^c^z 


X 


07  +  0ox^(0ox 

CP,  f 


07)  +  ^&(06)h2O 


^(07)h2O  • 

(22) 


2.4.  Electrochemical  model 


The  thermal  model  described  by  Eqs.  (12),  (13),  (18)-(22) 
needs  to  be  complemented  by  the  electrochemical  counterpart 
such  that  the  PEM  fuel  cell  power  generation  performance  can 
be  calculated.  The  resistances,  /3 (£?),  of  each  compartment  are 
evaluated  from  the  material  and  geometric  features  as 

Pi  =  ■  1/3  - - ,  j  =1,2.  6,  7  (23) 

Asvf  Oj0-4>j) 

Pj  =  ^~k - ,  j  =  3,  4,  5,  04  =  1  (24) 

AsVj  cfj(pj 

where  the  ionic  conductivity,  o  (£2_1  m_1),  of  Nafion  117  as 
a  function  of  temperature  is  given  by  the  following  empirical 
formula  [30] : 


o jiff)  =  exp 


7  =  3,  4,  5 


(0.5139A.J 


0.326), 


(25) 


The  conductivities  of  the  catalyst  layers  are  given  by  <73  03  and 
o 505,  according  to  Eqs.  (24)  and  (25),  which  agree  qualitatively 
with  previously  measured  catalyst  layers  ionic  conductivities 
[31],  i.e.,  the  ionic  conductivity  increases  with  increasing  Nafion 
content,  which  increases  as  0  increases  in  the  present  model.  The 
conductivities  of  the  diffusive  layers,  02  and  <76,  are  the  carbon- 
phase  conductivities  [32].  Finally,  the  conductivities  of  CV 1  and 
CV7,  o\  and  <r7,  are  the  electrical  conductivities  of  the  bipolar 
plate  material.  The  CV  1  and  CV7  void  fractions  in  Eq.  (23),  i.e., 
01  =  07  =  §ch Kj$t  +  §ch)>  are  computed  for  any  particular  fuel  cell 
internal  geometry  according  to  Fig.  1. 

In  Eq.  (25),  A  is  the  water  content  that  is  described  as  the  ratio 
of  the  number  of  water  molecules  to  the  number  of  charge  sites. 
In  a  polymer  membrane,  the  charge  sites  are  ions,  S03-H+.  It  is 
assumed  that  all  ingredients  of  the  catalyst  layer  are  evenly  dis¬ 
tributed,  and  that  the  liquid  water  product  is  evenly  distributed. 
Therefore,  the  water  content  in  ionomer  (within  the  catalyst 
layer)  is  assumed  constant:  this  corresponds  to  the  value  mea¬ 
sured  when  the  ionomer  is  in  contact  with  liquid  water,  or  to 
some  average  between  this  value  and  the  value  corresponding 
to  when  the  ionomer  is  in  contact  with  saturated  water  vapor 
[33].  Usually,  the  anode  water  content  is  different  from  that  in 
the  cathode;  therefore  for  assumed  values  of  Aa  (anode  water 
content,  or  A3)  and  Ac  (cathode  water  content,  or  A5),  and  by 
assuming  a  linear  variation  of  the  water  content  along  the  mem¬ 
brane  thickness,  the  average  water  content  in  the  membrane  is 
defined  as 

Aa  -j-  Xq 

A4  =  (26) 

Now,  we  can  calculate  the  molar  concentration  of  the  acid 
solution  [H(aq)+]  mentioned  in  CV3  energy  conservation  analy¬ 
sis  and  the  ionic  conductivity  in  Eq.  (25). 

The  polarization  curve  indicates  the  performance  of  a  fuel 
cell.  The  dimensionless  potential  is  defined  in  terms  of  a  given 
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reference  voltage  Vref ,  namely  V  =  V/  Vref  and  rj  =  r)/  Vref-  The 
dimensionless  actual  potential  V  is  an  accumulated  result  of 
dimensionless  anode  electrical  potential  Va,  dimensionless  cath¬ 
ode  electrical  potential  Vc,  and  the  dimensionless  ohmic  loss 
(rj ohm)  in  the  space  from  CV1  to  CV7,  i.e., 

V  =Va  +  Vc-  Oohm  (27) 


The  ohmic  loss  i)0hm  is  estimated  by 


Oohm  — 


(28) 


The  actual  electrical  potential  at  the  anode  is  ya  =  ye,a  — 
—  l^7d,aU  where  the  reversible  electrical  potential  at  the  anode 
ye,a  is  given  by  the  Nernst  equation  [29]: 


(29) 


where  Q3  =  [/7(aq)+]2/JPH2  and  V^a  =  2sG^/(—nF).  Mass  dif¬ 
fusion  and  charge  transfer  cause  the  potential  losses  at  the  anode. 
The  overpotential  due  to  charge  transfer  (r/a)  is  obtained  from 
the  Bulter-Volmer  equation  [34,35]  as  follows: 


exp 


f  (1  -Ua)0aF\ 

V  rt3  J 


UaOaF\ 

RT3  J 


(30) 


where  aa  is  the  anode  charge  transfer  coefficient,  and  /G, a  is 
the  anode  exchange  current  density  (a  function  of  catalyst  type, 
catalyst  layer  morphology,  temperature  and  pressure)  ,  which 
is  listed  in  Table  1.  The  potential  loss  due  to  mass  diffusion 
[34]  is 


RT3 

»?d,a  =  — tt  In 
nF 


^wet,30im,a 


(31) 


where  the  limiting  current  density  at  the  anode  ziim,a  occurs  at 
high  values  of  the  surface  overpotential,  when  the  gas  is  com¬ 
pletely  depleted  in  the  very  thin  active  catalyst  layer  fraction 
situated  at  the  interface  with  the  gas  diffuser,  i.e.,  T>2,out  =  0. 
Therefore,  /iim,a  is  obtained  from  Eq.  (15),  i.e., 


dim,  a 


PfD2nF 
47h2  L  2  Rf03  Tqq 


(32) 


The  methodology  in  estimating  the  anode  potential  is 
valid  in  building  the  cathode  potential  correlations.  Simi¬ 
larly,  the  actual  cathode  potential  is  Vc  =  ye,c  —  rjc  ~  \0d,c\ 
and  the  reversible  electrical  cathode  potential  is  Ve,c  = 

V'e’e  -  (RT5/nF)\n(Q5),  where  Q5  =  W Hiaql+ \2 pf)  1  and 
ye°c  =  A G°5/(—nF).  The  Butler- Volmer  equation  for  cal¬ 
culating  the  cathode  side  overpotential  ijc  is  7/^5, wet  = 
/0,c[exp((l  -  ac)y]cF/RT5)  -  exp(-acr]cF/RT5)].  The  cathode 

mass  diffusion  depleting  overpotential  is  7/d,c  =  RT3/nF  ln(l  — 

//Awet,5/iim,c)?  and  the  cathode  limiting  current  density  is 

dim,c  =  2/70x  D^nF/Mo2L^R 

ox  06Toc- 


2.5.  Pumping  power  and  fuel  cell  net  power  model 


The  dimensionless  pumping  power  Wp  is  required  to  supply 
the  fuel  cell  with  fuel  and  oxidant.  Therefore,  the  total  net  power 
output  (available  for  utilization)  of  the  fuel  cell  is 


Wnet  =  W  ~  %  (33) 

where  W  is  the  total  fuel  cell  electrical  power  output  and 
W  =  V L  It  is  necessary  to  evaluate  the  pressure  drops  in  the 
reactant  supply  channel  before  estimating  the  pumping  power 
consumption  in  the  fuel  cell.  Assuming  that  the  channel  is 
straight  and  sufficiently  slender,  the  pressure  drops  are  expressed 
by 


Q  j  Rk  Toe 


(34) 


where 7  =  1,7  and  k  =f  ox,  respectively.  The  friction  factor  (//)  in 
the  gas  channels  is  estimated  by  experimental  correlations  from 
the  literature.  For  the  laminar  regime  (Re h  <  2300)  we  used  the 
correlation  [33]: 


fj  RehJ  =  24(1  -  1.3553 Sj  +  1.9467^  -  1.70125] 

+0.95645]  -  0.25375])  (35) 

where  8j  =  Lc\l/Lj,  for  Tch  <  Lj  and  8j  =  Lj/Lc h,  for  Lc\l>Lj; 
l\i  =  2LchLj/(Lcb  +  Lj),  Re\j  =  u;7)h.;/>///(  /  and  7=1,7.  The  cor- 
relation  used  for  the  turbulent  regime  is  the  same  as  the  CV 1 
energy  conservation  analysis. 

Then,  the  dimensionless  pumping  power  Wp  is  obtained  as 
follows: 

Wp  =  VrfSf^-APi  +  ^ox S0x  -f  A  Pi  (36) 

P 1  Pi 

where 


fh  ref  Too  R  j 

yref  7ref 


(37) 


3.  Shape  optimization  and  results 


The  thermal  system  reduces  to  seven  nonlinear  ordinary 
differential  equations,  Eqs.  (12),  (13),  (18)-(22),  in  which 
the  unknowns  are  the  temperatures  of  the  seven  control  vol¬ 
umes.  A  Fortran  computational  code  was  written  to  obtain  both 
steady  state  and  transient  solutions.  The  transient  solutions  were 
obtained  by  solving  the  system  with  a  Runge-Kutta  fourth/fifth 
order  method  [36]  from  a  given  initial  condition.  For  obtaining 
steady  state  results,  a  system  of  seven  nonlinear  algebraic  equa¬ 
tions  was  generated  by  setting  the  time  derivatives  equal  to  0. 
The  system  was  then  solved  with  a  quasi-Newton  method  [36]. 
The  Newton  iterative  process  was  set  to  achieve  convergence 
when  the  Euclidean  norm  of  the  residual  of  the  system  was  less 
than  10-6.  Because  of  this,  the  computational  time  required  for 
one  geometric  structure  study  case  was  short. 

In  order  to  make  the  PEM  fuel  cell  design  close  to  the  vehicle 
application,  in  this  study  we  conducted  a  multi-objective  opti¬ 
mization  procedure  to  seek  the  optimal  geometric  dimension  or 


364 


S.  Chen  et  al.  /  Journal  of  Power  Sources  162  (2006)  356-368 


limits  with  the  goal  of  immediate  response  to  step  current  load 
and  as  high  as  possible  steady- state  net  power  output.  Table  1 
shows  the  properties  and  constant  values  used  to  produce  the 
numerical  results  of  the  present  study.  The  fuel  and  oxidant 
inlet  temperature  were  assumed  at  298.15  K.  The  average  fuel 
and  oxidant  channel  pressures  were  assumed  as  1  and  1.2  atm, 
respectively.  Fig.  1  illustrates  the  main  geometric  features  of  the 
single  PEM  fuel  considered  in  this  study:  the  external  dimen¬ 
sions  (Lx,  Ly ,  and  Lz ),  and  the  individual  layer  thicknesses  (L/). 
The  geometric  constraints  of  the  original  study  case  of  this  paper 
are  as  follows,  VT  =  2.25,  fy  =  §z  =  160£*,  and  (£i,  §3,  §4,  §5, 

§6,  §7)  =  (0.05,  0.36,  0.01,  0.06,  0.01,  0.36,  0.05)^. 

The  PEM  fuel  cell  polarization  curve  indicates  that  there 
exists  an  optimal  working  current  at  which  the  fuel  cell  has 
maximal  net  power  output.  In  this  study,  we  focus  on  the  fuel 
cell  steady-state  performance  (net  power  output  Wnet,s)  and  the 
transient  behavior  (start-up  time  tsu )  characteristics  under  such 
optimal  working  current  condition.  We  also  explore  the  con¬ 
sidered  geometric  dimension  effects  on  those  two  performance 
indicators.  Fig.  2a  and  b  illustrate  a  typical  transient  behavior  of  a 
PEM  fuel  cell  in  its  start-up  process.  Fig.  2a  shows  the  electrode 
potentials  and  the  power  output  compartments  in  the  transient 
process  responding  to  a  step-changed  current  load.  Fig.  2b  shows 
the  dimensionless  temperature  transient  process  after  a  step- 
changed  current  load  applied  to  the  fuel  cell.  Considering  the 
power  support  purpose  of  the  PEM  fuel  cell,  the  transient  behav¬ 
ior  indicator,  the  start-up  time,  is  defined  as  the  time  when  the 
net  power  output  arrived  at  98%  of  its  steady-state  power  output. 
This  measurement  takes  into  account  two  phenomena  of  PEM 
fuel  cell  transient  behavior  shown  in  Fig.  2a  and  b:  (i)  electri¬ 
cal  properties  (e.g.  voltages,  net  power  output)  respond  more 
quickly  than  the  thermal  properties  (e.g.  temperature)  -  there 
are  power  and  voltage  jumps  in  the  beginning  of  the  transient 
process,  while  the  temperatures  of  each  component  of  fuel  cell 
increase  slowly  and  continuously  because  of  high  thermal  iner¬ 
tia  of  porous  media  and  the  reactants  -  this  phenomenon  is  quite 
similar  to  the  conventional  internal  combustion  engine  that  sup¬ 
plies  the  power  requirement  almost  immediately,  but  arrives  to 
its  thermal  steady  state  quite  slowly,  and  (ii)  after  the  net  power 
output  reaches  98%  of  its  steady-state  power  level,  the  net  power 
output  curve’s  slope  is  close  to  0,  while  the  fuel  cell  temperatures 
are  still  increasing — in  other  words,  after  a  certain  period  of  the 
transient  process,  the  thermal  mechanisms  (e.g.  temperatures) 
do  not  affect  significantly  the  PEM  fuel  cell  power  generation. 
Therefore,  it  is  questionable  to  consider  the  transient  period  as 
the  start-up  time  based  on  temperature  stabilization. 

Fig.  3  demonstrates  that  the  electrode  reaction  layer  thick¬ 
nesses  (§3  and  §5)  are  crucial  to  PEM  fuel  cell  transient  pro¬ 
cess  which  can  be  categorized  into  three  modes,  i.e.  “under¬ 
powered”,  “exact-powered”,  and  “over-powered”.  When  the 
electrode  reaction  thickness  is  small,  say  smaller  than  0.003^, 
the  transient  process  is  “over-powered”,  while  the  net  power 
output  jumps  higher  than  the  steady-state  net  power  output  and 
then  decreases  to  the  steady-state  net  power  output  level.  When 
the  electrode  reaction  layer  thickness  is  increased  to  0.007^, 
the  net  power  output  reaches  the  steady-state  level  immediately, 
named  as  “exact-powered”  transient  process.  When  the  electrode 
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Fig.  2.  Polarization  and  power  curves  with  respect  to  time  (a),  and  net  power 
and  internal  temperature  vs.  time  (b). 


reaction  layer  thickness  exceeds  this  critical  value,  the  tran¬ 
sient  process  will  be  “under-powered”  and  the  net  power  output 
increases  to  the  steady- state  level  continuously  after  the  begin¬ 
ning  jump.  Fig.  3  was  constructed  in  a  way  to  show  clearly  the 
three  transient  operation  modes,  but  it  should  be  noted  that  Wnet,s 
is  strongly  dependent  on  the  reaction  layer  thickness,  which  is 
the  subject  of  discussion  in  Fig.  4. 

Fig.  4  shows  the  comparison  of  the  PEM  fuel  cell  steady-state 
performance,  the  steady-state  net  power  output,  and  the  tran¬ 
sient  behavior,  i.e.,  the  start-up  time.  The  PEM  fuel  cell  shows 
maximum  steady- state  net  power  output  and  short  start-up  time 
of  “over-powered”  process  at  §3  =  §5  =0.003^.  It  is  the  multi¬ 
objective  optimization  design  of  the  electrode  reaction  layer  for 
the  PEM  fuel  cell  and  the  first  level  of  optimization  in  this  study, 
named  one-way  optimization.  In  the  “under-powered”  range,  §3, 
§5  >0.007^,  the  steady-state  net  power  output  decreases  dra¬ 
matically,  and  the  start-up  time  increases  considerably,  therefore 
affecting  negatively  the  two  objective  functions. 
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Fig.  3.  The  normalized  net  power  output  of  PEMFC  transient  response. 

A  fixed  internal  geometric  parameter  is  the  thickness  (§4)  of 
the  polymer  electrolyte  membrane.  The  optimization  of  the  elec¬ 
trode  reaction  layers  thickness  was  conducted  for  several  values 
of  §4,  which  produces  the  results  shown  in  Fig.  5.  It  is  observed 
that  the  transient  process  is  in  the  “over-powered”  or  “exact- 
powered”  modes,  and  the  start-up  time  is  close  to  0,  when  ^4/^ 
is  smaller  than  0.08.  When  the  membrane  thickness  exceeds  such 
critical  ratio,  the  transient  process  deteriorates  dramatically,  and 
the  start-up  time  increases  abruptly.  In  the  entire  range  of  the 
membrane  thickness  studied,  the  steady-state  net  power  output 
decreases  monotonically  as  the  membrane  thickness  increases 
because  of  the  increment  of  proton  transport  resistance  in  the 
membrane  and  the  drop  of  the  membrane  ionic  conductivity. 
The  start-up  time  curve  is  unimodal  and  approaches  0  when 
H\!%x  <  0.08  due  to  the  small  electrical  inertia  and  the  low  elec¬ 
trical  resistance.  When  §4/^  >0.08,  the  start-up  time  increases 
drastically  and  reaches  a  maximum  value  at  §4 /§*  =  0.2  because 


Fig.  5.  The  effect  of  membrane  thickness  on  the  one-way  maximized  net  power. 

as  the  membrane  thickness  increases,  the  fuel  cell  thermal  and 
electrical  inertias  increase  as  well.  When  §4 /§*  >0.2,  the  optimal 
working  current  decreases  as  the  membrane  thickness  increases, 
which  in  turn  reduces  the  start-up  time.  Fig.  5  also  indicates  that 
the  membrane  should  be  as  thin  as  possible,  down  to  its  minimum 
dielectric  limit,  which  is  also  restricted  by  the  manufacturing 
technology.  So,  although  Fig.  4  shows  that  rsu  increases  mono¬ 
tonically  as  the  reaction  layers  thickness  increase,  it  also  shows 
the  one-way  optimized  structure  leads  to  short  start-up  time 
for  a  given  membrane  thickness.  As  the  membrane  thickness 
increases,  even  with  the  one-way  optimized  external  structure, 
rsu  increases,  which  is  shown  in  Fig.  5.  The  useful  design  infor¬ 
mation  is  that  ^4/^x  <  0.08  for  short  start-up  time. 

The  previous  internal  structure  optimization  process  is 
repeated  at  different  external  aspect  ratios  to  produce  the  results 
shown  in  Fig.  6.  We  find  that  the  start-up  time  curve  departs 
from  the  values  close  to  0  when  the  external  aspect  ratios, 
^y/^x  and  are  smaller  than  120,  and  ^y/^x  =  HJHx  >  0.001. 
The  steady- state  net  power  output  decreases  with  the  incre¬ 
ment  of  external  aspect  ratios  within  the  considered  range 
Wx  =  HJHx  >  100).  Other  external  aspect  ratio  cases  are  also 
studied  and  showed  in  Fig.  7.  The  steady- state  net  power 
output  reaches  a  maximum  when  the  external  aspect  ratio  is 
approximately  90,  ^y/^x  =  HJHx  >  90.  However,  the  PEM  fuel 
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Fig.  4.  The  one-way  net  power  maximization  for  optimal  internal  structure. 


Fig.  6.  The  effect  of  external  shape  on  net  power  and  start-up  time. 
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Fig.  7.  The  two-way  net  power  maximization  for  optimal  internal  structure  and 
external  shape. 


cell  responds  to  the  external  step  current  loads  rather  slowly. 
The  start-up  time  exceeds  10  s  for  the  parameters  of  Table  1. 
In  order  to  guarantee  a  short  response  (or  start-up)  time,  one 
should  increase  the  external  aspect  ratio  at  the  expense  of 
reducing  the  maximized  steady-state  net  power  output.  A  value 
of  ^y/^x  =  ^z/^x>  120  leads  to  tsu  =  30  (or  fsu^0.1s),  with  a 
decrease  of  only  9%  in  the  maximized  net  power  output. 

Several  different  total  volumes  over  the  range  1.0-10.0,  for 
the  stoichiometric  ratios  =  §7  =  2,  are  studied  in  Figs.  8  and  9. 
From  Fig.  8,  it  is  observed  that  there  exists  a  critical  external 
aspect  ratio  for  each  volume.  The  start-up  time  increases 
drastically  when  the  external  aspect  ratio  drops  below  those 
critical  values.  Fig.  9  shows  the  two-way  multi-objective  optima 
in  a  log-log  graph,  adopting  the  critical  values  for  the  external 
aspect  ratio.  The  optimal  external  aspect  ratios  are  propor- 
tional  to  Vf2497,  i.e.,  (£,/&  =  fe/^)opt  =  94.3843 Vf2497, 
and  the  two-way  maximum  steady- state  net  power  output 
associated  with  immediate  response  to  the  external  load  is 
proportional  to  V^-6321,  i.e.,  lTnet,s,mm  =  155.2745 V^-6321. 
The  effects  of  different  stoichiometric  ratios  ranging  from 
1.0  to  10.0,  on  the  maximum  steady-state  net  power  output 
and  on  the  optimal  external  aspect  ratio  associated  with  the 
two-way  multi-objective  optimization  along  the  several  total 
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Fig.  8.  The  external  shape  effect  on  start-up  time. 


Fig.  9.  The  behavior  of  the  optimal  shape  and  the  two-way  maximized  net  power 
output  with  respect  to  total  volume. 


(b)  V 


Fig.  10.  The  effect  of  stoichiometric  ratio  and  total  volume  on  the  behavior  of 
the  two-way  maximized  net  power  (a)  and  optimal  shape  (b). 

volumes  over  the  range  1.0-10.0  are  shown  in  Fig.  10a  and 
b,  respectively.  The  two-way  maximized  net  power  output 
fitting  equations  are  Wnet  s  mm  =  297.2350Vif6432  when 
ft  =  ?7  =  1,  Wnet,s,mm  =  155.2744  V°'6321  when  f1=f1  = 
2,  and  Wnet,  s,  mm  =  39.2103  Vj 6239  when  £i  =  =  10. 
The  optimal  aspect  ratio  fitting  equations  are 
%!Hx  =  fe/^)opt  =  145.2446 Vf2616  when  O  =  Kl  =  1. 
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%/Hx  =  fe/^)opt  =  94.3843 yo-2497  when  ?i=?7  =  2,  and 

(&/&  =  fe/^)opt  =  40.4390  Vf2409  when  ?i  =  ?7  =  10.  The 
optimal  external  aspect  ratio  is  approximately  proportion  to 

VT  •  The  two-way  maximized  net  power  output  is  approxi¬ 
mately  proportional  to  V^  63.  This  power-volume  correlation 
is  also  found  in  nature.  According  to  Randall  et  al.  [37],  the 
metabolic  rate  within  a  single  species  is  around  0.63  instead 
of  2/3,  because  of  the  thickening  of  bones  and  other  body 
structures  in  larger  animals. 

4.  Conclusions 

In  this  paper,  we  showed  a  multi-objective  optimization  of 
a  single  PEM  fuel  cell  concentrating  on  vehicle  applications, 
which  require  the  PEM  fuel  cell  to  respond  to  the  external  load 
quickly  and  to  provide  as  much  power  as  possible.  The  internal 
and  external  shape  effects  on  the  PEM  fuel  cell  steady  state  and 
transient  response  were  investigated.  The  trade-off  between  the 
steady  state  net  power  output  and  the  start-up  time  was  optimally 
balanced  in  the  multi-objective  optimization.  The  geometric 
degrees  of  freedom  were  identified  physically  to  be  optimized, 
i.e.,  the  electrodes  reaction  layers  thickness,  L3  and  L5,  and  the 
external  aspect  ratios,  LyILx  and  Lz/Lx.  The  optima  found  are 
sharp  and  therefore,  important  to  be  identified  in  actual  PEM 
fuel  cell  design. 

The  key  conclusion  is  that  we  showed  that  the  trade-offs  exist, 
and  that  from  them  results  the  internal  and  external  structure  of 
a  single  PEM  fuel  cell,  i.e.,  constructal  design  [20,21,22,38]. 
The  maximum  net  power  output  is  proportional  to  V^-63.  This 
exponent  is  the  metabolic  rate  within  a  single  species  in  nature. 

The  optimal  external  aspect  ratio  is  proportional  to  VjA.  These 
two  relations  are  crucial  to  the  auxiliary  power  units  systems 
(APUS)  design.  In  principle,  the  multi- objective  optimization 
illustrated  in  this  study  can  be  extended  on  a  hierarchical  ladder 
to  large  and  more  complex  PEM  fuel  cell  systems,  to  explore 
multi-scale  packings  that  use  the  available  volume  for  maximum 
power  generation  and  quick  response. 
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